In [1]:
import sys
!{sys.executable} -m pip install -r requirements.txt
Collecting cvxpy==1.0.3 (from -r requirements.txt (line 1))
  Downloading https://files.pythonhosted.org/packages/a1/59/2613468ffbbe3a818934d06b81b9f4877fe054afbf4f99d2f43f398a0b34/cvxpy-1.0.3.tar.gz (880kB)
    100% |████████████████████████████████| 880kB 13.9MB/s ta 0:00:01    73% |███████████████████████▌        | 645kB 25.5MB/s eta 0:00:01
Collecting numpy==1.16.0 (from -r requirements.txt (line 2))
  Downloading https://files.pythonhosted.org/packages/7b/74/54c5f9bb9bd4dae27a61ec1b39076a39d359b3fb7ba15da79ef23858a9d8/numpy-1.16.0-cp36-cp36m-manylinux1_x86_64.whl (17.3MB)
    100% |████████████████████████████████| 17.3MB 2.2MB/s eta 0:00:01   20% |██████▊                         | 3.6MB 28.4MB/s eta 0:00:01    28% |█████████▏                      | 5.0MB 22.9MB/s eta 0:00:01    35% |███████████▍                    | 6.2MB 28.3MB/s eta 0:00:01    43% |█████████████▉                  | 7.5MB 28.8MB/s eta 0:00:01    95% |██████████████████████████████▋ | 16.5MB 26.7MB/s eta 0:00:01
Collecting pandas==0.21.1 (from -r requirements.txt (line 3))
  Downloading https://files.pythonhosted.org/packages/3a/e1/6c514df670b887c77838ab856f57783c07e8760f2e3d5939203a39735e0e/pandas-0.21.1-cp36-cp36m-manylinux1_x86_64.whl (26.2MB)
    100% |████████████████████████████████| 26.2MB 1.4MB/s eta 0:00:01  0% |▎                               | 245kB 23.3MB/s eta 0:00:02    39% |████████████▊                   | 10.4MB 25.9MB/s eta 0:00:01    58% |██████████████████▋             | 15.3MB 25.7MB/s eta 0:00:01    71% |███████████████████████         | 18.8MB 24.8MB/s eta 0:00:01    88% |████████████████████████████▌   | 23.3MB 22.3MB/s eta 0:00:01    93% |█████████████████████████████▊  | 24.4MB 24.7MB/s eta 0:00:01    96% |███████████████████████████████ | 25.4MB 24.0MB/s eta 0:00:01
Collecting plotly==2.2.3 (from -r requirements.txt (line 4))
  Downloading https://files.pythonhosted.org/packages/99/a6/8214b6564bf4ace9bec8a26e7f89832792be582c042c47c912d3201328a0/plotly-2.2.3.tar.gz (1.1MB)
    100% |████████████████████████████████| 1.1MB 15.0MB/s ta 0:00:01
Collecting scipy==1.0.0 (from -r requirements.txt (line 5))
  Downloading https://files.pythonhosted.org/packages/d8/5e/caa01ba7be11600b6a9d39265440d7b3be3d69206da887c42bef049521f2/scipy-1.0.0-cp36-cp36m-manylinux1_x86_64.whl (50.0MB)
    100% |████████████████████████████████| 50.0MB 787kB/s eta 0:00:01  2% |▋                               | 1.0MB 22.6MB/s eta 0:00:03    6% |██                              | 3.2MB 23.5MB/s eta 0:00:02    8% |██▉                             | 4.4MB 23.6MB/s eta 0:00:02    12% |████▏                           | 6.5MB 23.6MB/s eta 0:00:02    17% |█████▌                          | 8.6MB 23.9MB/s eta 0:00:02    19% |██████▏                         | 9.7MB 21.2MB/s eta 0:00:02    21% |███████                         | 10.9MB 22.4MB/s eta 0:00:02    23% |███████▋                        | 12.0MB 22.2MB/s eta 0:00:02    25% |████████▎                       | 13.0MB 23.5MB/s eta 0:00:02    32% |██████████▍                     | 16.3MB 22.6MB/s eta 0:00:02    34% |███████████                     | 17.3MB 26.3MB/s eta 0:00:02    36% |███████████▉                    | 18.4MB 22.3MB/s eta 0:00:02    38% |████████████▍                   | 19.4MB 22.1MB/s eta 0:00:02    41% |█████████████▏                  | 20.5MB 23.1MB/s eta 0:00:02    43% |█████████████▉                  | 21.6MB 23.0MB/s eta 0:00:02    45% |██████████████▋                 | 22.8MB 24.8MB/s eta 0:00:02    49% |████████████████                | 25.0MB 25.2MB/s eta 0:00:01    52% |████████████████▊               | 26.1MB 18.2MB/s eta 0:00:02    57% |██████████████████▎             | 28.5MB 26.5MB/s eta 0:00:01    61% |███████████████████▋            | 30.7MB 27.7MB/s eta 0:00:01    63% |████████████████████▎           | 31.8MB 22.1MB/s eta 0:00:01    70% |██████████████████████▍         | 35.1MB 23.5MB/s eta 0:00:01    74% |███████████████████████▉        | 37.2MB 27.1MB/s eta 0:00:01    76% |████████████████████████▋       | 38.5MB 23.8MB/s eta 0:00:01    80% |█████████████████████████▉      | 40.4MB 20.3MB/s eta 0:00:01    82% |██████████████████████████▌     | 41.5MB 24.5MB/s eta 0:00:01    91% |█████████████████████████████▏  | 45.5MB 19.4MB/s eta 0:00:01    95% |██████████████████████████████▍ | 47.6MB 29.5MB/s eta 0:00:01    99% |███████████████████████████████▊| 49.6MB 28.0MB/s eta 0:00:01
Collecting tensorflow-tensorboard==0.1.8 (from -r requirements.txt (line 6))
  Downloading https://files.pythonhosted.org/packages/93/31/bb4111c3141d22bd7b2b553a26aa0c1863c86cb723919e5bd7847b3de4fc/tensorflow_tensorboard-0.1.8-py3-none-any.whl (1.6MB)
    100% |████████████████████████████████| 1.6MB 14.3MB/s ta 0:00:01    18% |██████                          | 307kB 28.3MB/s eta 0:00:01
Collecting osqp (from cvxpy==1.0.3->-r requirements.txt (line 1))
  Downloading https://files.pythonhosted.org/packages/de/90/aee4a9a72bd8643c0dd2cbc1735c5f4fa9fc6c0131935b60cb167eaeb18b/osqp-0.6.2.post5-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl (232kB)
    100% |████████████████████████████████| 235kB 10.3MB/s ta 0:00:01    52% |█████████████████               | 122kB 12.1MB/s eta 0:00:01
Collecting ecos>=2 (from cvxpy==1.0.3->-r requirements.txt (line 1))
  Downloading https://files.pythonhosted.org/packages/08/0c/a9195eb04d6b58572a4379d98661d138066da97fb8dbbdc4e104452e8377/ecos-2.0.10.tar.gz (135kB)
    100% |████████████████████████████████| 143kB 14.8MB/s ta 0:00:01
  Installing build dependencies ... done
Collecting scs>=1.1.3 (from cvxpy==1.0.3->-r requirements.txt (line 1))
  Downloading https://files.pythonhosted.org/packages/f0/a4/e3cbddd60818e7e900fd0a8e0f0555692c353be7aa8ca1a8299a05c20699/scs-3.2.0.tar.gz (657kB)
    100% |████████████████████████████████| 665kB 17.4MB/s ta 0:00:01
  Installing build dependencies ... done
Collecting multiprocess (from cvxpy==1.0.3->-r requirements.txt (line 1))
  Downloading https://files.pythonhosted.org/packages/df/d0/c3011a5bb77f307e68682a5046cce1a2c6591267bf24b5bf3fc4130bb39d/multiprocess-0.70.12.2-py36-none-any.whl (106kB)
    100% |████████████████████████████████| 112kB 13.0MB/s ta 0:00:01
Requirement already satisfied: fastcache in /opt/conda/lib/python3.6/site-packages (from cvxpy==1.0.3->-r requirements.txt (line 1)) (1.0.2)
Requirement already satisfied: six in /opt/conda/lib/python3.6/site-packages (from cvxpy==1.0.3->-r requirements.txt (line 1)) (1.11.0)
Requirement already satisfied: toolz in /opt/conda/lib/python3.6/site-packages (from cvxpy==1.0.3->-r requirements.txt (line 1)) (0.8.2)
Requirement already satisfied: python-dateutil>=2 in /opt/conda/lib/python3.6/site-packages (from pandas==0.21.1->-r requirements.txt (line 3)) (2.6.1)
Requirement already satisfied: pytz>=2011k in /opt/conda/lib/python3.6/site-packages (from pandas==0.21.1->-r requirements.txt (line 3)) (2017.3)
Requirement already satisfied: decorator>=4.0.6 in /opt/conda/lib/python3.6/site-packages (from plotly==2.2.3->-r requirements.txt (line 4)) (4.0.11)
Requirement already satisfied: nbformat>=4.2 in /opt/conda/lib/python3.6/site-packages (from plotly==2.2.3->-r requirements.txt (line 4)) (4.4.0)
Requirement already satisfied: requests in /opt/conda/lib/python3.6/site-packages (from plotly==2.2.3->-r requirements.txt (line 4)) (2.18.4)
Requirement already satisfied: bleach==1.5.0 in /opt/conda/lib/python3.6/site-packages (from tensorflow-tensorboard==0.1.8->-r requirements.txt (line 6)) (1.5.0)
Requirement already satisfied: werkzeug>=0.11.10 in /opt/conda/lib/python3.6/site-packages (from tensorflow-tensorboard==0.1.8->-r requirements.txt (line 6)) (0.14.1)
Requirement already satisfied: wheel>=0.26 in /opt/conda/lib/python3.6/site-packages (from tensorflow-tensorboard==0.1.8->-r requirements.txt (line 6)) (0.30.0)
Requirement already satisfied: protobuf>=3.2.0 in /opt/conda/lib/python3.6/site-packages (from tensorflow-tensorboard==0.1.8->-r requirements.txt (line 6)) (3.5.1)
Requirement already satisfied: html5lib==0.9999999 in /opt/conda/lib/python3.6/site-packages (from tensorflow-tensorboard==0.1.8->-r requirements.txt (line 6)) (0.9999999)
Requirement already satisfied: markdown>=2.6.8 in /opt/conda/lib/python3.6/site-packages (from tensorflow-tensorboard==0.1.8->-r requirements.txt (line 6)) (2.6.9)
Collecting qdldl (from osqp->cvxpy==1.0.3->-r requirements.txt (line 1))
  Downloading https://files.pythonhosted.org/packages/08/8c/35bf914d768e29954617113597ed313340da730f4a8a58374a0dc70295c5/qdldl-0.1.5.post2.tar.gz (69kB)
    100% |████████████████████████████████| 71kB 10.7MB/s ta 0:00:01
Collecting dill>=0.3.4 (from multiprocess->cvxpy==1.0.3->-r requirements.txt (line 1))
  Downloading https://files.pythonhosted.org/packages/b6/c3/973676ceb86b60835bb3978c6db67a5dc06be6cfdbd14ef0f5a13e3fc9fd/dill-0.3.4-py2.py3-none-any.whl (86kB)
    100% |████████████████████████████████| 92kB 11.2MB/s ta 0:00:01
Requirement already satisfied: ipython-genutils in /opt/conda/lib/python3.6/site-packages (from nbformat>=4.2->plotly==2.2.3->-r requirements.txt (line 4)) (0.2.0)
Requirement already satisfied: jupyter-core in /opt/conda/lib/python3.6/site-packages (from nbformat>=4.2->plotly==2.2.3->-r requirements.txt (line 4)) (4.4.0)
Requirement already satisfied: jsonschema!=2.5.0,>=2.4 in /opt/conda/lib/python3.6/site-packages (from nbformat>=4.2->plotly==2.2.3->-r requirements.txt (line 4)) (2.6.0)
Requirement already satisfied: traitlets>=4.1 in /opt/conda/lib/python3.6/site-packages (from nbformat>=4.2->plotly==2.2.3->-r requirements.txt (line 4)) (4.3.2)
Requirement already satisfied: chardet<3.1.0,>=3.0.2 in /opt/conda/lib/python3.6/site-packages (from requests->plotly==2.2.3->-r requirements.txt (line 4)) (3.0.4)
Requirement already satisfied: idna<2.7,>=2.5 in /opt/conda/lib/python3.6/site-packages (from requests->plotly==2.2.3->-r requirements.txt (line 4)) (2.6)
Requirement already satisfied: urllib3<1.23,>=1.21.1 in /opt/conda/lib/python3.6/site-packages (from requests->plotly==2.2.3->-r requirements.txt (line 4)) (1.22)
Requirement already satisfied: certifi>=2017.4.17 in /opt/conda/lib/python3.6/site-packages (from requests->plotly==2.2.3->-r requirements.txt (line 4)) (2019.11.28)
Requirement already satisfied: setuptools in /opt/conda/lib/python3.6/site-packages (from protobuf>=3.2.0->tensorflow-tensorboard==0.1.8->-r requirements.txt (line 6)) (38.4.0)
Building wheels for collected packages: cvxpy, plotly, ecos, scs, qdldl
  Running setup.py bdist_wheel for cvxpy ... done
  Stored in directory: /root/.cache/pip/wheels/2b/60/0b/0c2596528665e21d698d6f84a3406c52044c7b4ca6ac737cf3
  Running setup.py bdist_wheel for plotly ... done
  Stored in directory: /root/.cache/pip/wheels/98/54/81/dd92d5b0858fac680cd7bdb8800eb26c001dd9f5dc8b1bc0ba
  Running setup.py bdist_wheel for ecos ... done
  Stored in directory: /root/.cache/pip/wheels/34/83/48/2192df21a678bf0fde91e59ec6ff78525bb3d99919397d6a8c
  Running setup.py bdist_wheel for scs ... done
  Stored in directory: /root/.cache/pip/wheels/3e/22/88/9a0da5b5c39fb5a10ad2e7c6192a268487d117770e24233980
  Running setup.py bdist_wheel for qdldl ... done
  Stored in directory: /root/.cache/pip/wheels/ad/28/83/0674c1f775e91a57943cf178dedf3a40cf48d9bf63209bd76c
Successfully built cvxpy plotly ecos scs qdldl
Installing collected packages: numpy, scipy, qdldl, osqp, ecos, scs, dill, multiprocess, cvxpy, pandas, plotly, tensorflow-tensorboard
  Found existing installation: numpy 1.12.1
    Uninstalling numpy-1.12.1:
      Successfully uninstalled numpy-1.12.1
  Found existing installation: scipy 1.2.1
    Uninstalling scipy-1.2.1:
      Successfully uninstalled scipy-1.2.1
  Found existing installation: dill 0.2.7.1
    Uninstalling dill-0.2.7.1:
      Successfully uninstalled dill-0.2.7.1
  Found existing installation: pandas 0.23.3
    Uninstalling pandas-0.23.3:
      Successfully uninstalled pandas-0.23.3
  Found existing installation: plotly 2.0.15
    Uninstalling plotly-2.0.15:
      Successfully uninstalled plotly-2.0.15
Successfully installed cvxpy-1.0.3 dill-0.3.4 ecos-2.0.10 multiprocess-0.70.12.2 numpy-1.16.0 osqp-0.6.2.post5 pandas-0.21.1 plotly-2.2.3 qdldl-0.1.5.post2 scipy-1.0.0 scs-3.2.0 tensorflow-tensorboard-0.1.8
In [2]:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
import cvxpy as cvx
import plotly as py
import plotly.graph_objs as go
import helper

py.offline.init_notebook_mode(connected=True)
%matplotlib inline
plt.style.use('ggplot')

Optimization with an Alpha Model and a Risk Model

In this exercise, we will construct an optimization problem using a stock universe consisting of 3 stocks. This problem will inherently be somewhat artificial, but using a stock universe of 3 stocks will allow us to create plots in 3 dimensions. These are useful for illustrating exactly what is going on in each step and will help develop your intuition for the problem. In this exercise, we will do the following steps:

  1. Create a 1-year-momentum-based alpha vector from some stock price data.
  2. Create a risk model using PCA.
  3. Construct the optimization problem using the alpha vector and risk model, and apply a set of standard constraints.
  4. Solve the optimization problem.

Load the data

Load the data from the file stock_prices_advanced_optimization.csv.

In [3]:
prices = pd.read_csv('stock_prices_advanced_optimization.csv', parse_dates=['date'], index_col=0)

The data are the price trends for 3 stocks, Stock A, Stock B and Stock C, for 4 years, from 2013 to 2017.

In [4]:
prices.head()
Out[4]:
A B C
date
2013-07-02 100.000000 102.000000 104.000000
2013-07-03 99.872139 102.298629 104.821312
2013-07-05 99.109975 102.733517 104.256297
2013-07-08 99.971278 102.877047 103.631199
2013-07-09 98.866192 102.568696 106.375325
In [5]:
prices.plot()
Out[5]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f1933c3fe48>

Calculate the returns of these price data.

In [6]:
rets = prices.pct_change()[1:].fillna(0)

Create an alpha vector

In [7]:
from scipy.stats import zscore

We'll create an alpha vector based on 1-year momentum.

In [8]:
# 1-yr momentum alpha

def log_returns(series):
    return np.log(series[-1])-np.log(series[0])
    
alpha = prices.rolling(window=252).apply(log_returns).rank(axis='columns').apply(zscore, axis='columns')

Now we'll take the row of most recent alpha values to use as our alpha vector. There should be a value in the vector for each stock.

In [9]:
# Take most recent set of values
alpha_vector = alpha.iloc[-1,:]
print(alpha_vector)
A    0.000000
B    1.224745
C   -1.224745
Name: 2017-06-30 00:00:00, dtype: float64

For this problem, we'll use the optimization objective $-\boldsymbol{\alpha}^T \mathbf{x}$. Remember, we are trying to minimize this function (to maximize alpha). Let's plot $-\boldsymbol{\alpha}^T \mathbf{x}$ as a function of the components of $\mathbf{x}$ so that we can get a better sense of the "shape" of the function.

In [10]:
n = 10
x = y = z = np.linspace(-2,2,n)

xv, yv, zv = np.meshgrid(x, y, z, indexing='ij')
obj_val = np.full(xv.shape, np.nan)

for i in range(n):
    for j in range(n):
        for k in range(n):
            obj_val[i,j,k] = -alpha_vector[0]*xv[i,j,k]-alpha_vector[1]*yv[i,j,k]-alpha_vector[2]*zv[i,j,k]
In [11]:
hover_text = helper.generate_hover_text(xv, yv, zv, 'Weight on Stock A', 'Weight on Stock B', 'Weight on Stock C', fcn_values=obj_val, fcn_label='Objective Function')

trace1 = go.Scatter3d(
    x=xv.flatten(),
    y=yv.flatten(),
    z=zv.flatten(),
    text = hover_text.flatten(),
    mode='markers',
    marker=dict(
        size=4,
        color=obj_val.flatten(),     # set color to an array/list of desired values
        colorscale='Viridis',   # choose a colorscale
        colorbar=dict(
                title='Objective Function'
            ),
        opacity=0.4,
        showscale=True
    ),
    hoverinfo = 'text'
)

data = [trace1]
layout = helper.create_standard_layout('Alpha Vector Optimization Objective Function', 'Weight on')
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig)

In this plot, the value of the objective function is represented by the color at each point in weight space. Notice that the function increases with weight on Stock C, decreases with weight on Stock B, and doesn't depend on the weight on Stock A.

What happens if we try to run the optimization now, by trying to minimize this objective function, under no constraints?

In [12]:
def find_optimal_holdings(alpha_vector):
    """
    Use cvxpy to construct and solve an optimization problem. Use -alpha*x as the objective, and use no constraints.

    Parameters
    ----------
    alpha_vector : DataFrame
        The 3-stock alpha vector calculated above.

    Returns
    -------
    optimal_weights : DataFrame
        A DataFrame containing the optimal weights calculated by the optimization algorithm.
    """
    #TODO: Implement function
    m = len(alpha_vector)
    x = cvx.Variable(m)
    
    obj = cvx.Minimize(-alpha_vector.values.flatten()*x)
    constraints=[]
    
    prob = cvx.Problem(obj, constraints)
    prob.solve(verbose=True, max_iter=500)
    
    optimal_weights = np.asarray(x.value).flatten()
    
    return pd.DataFrame(data=optimal_weights)
In [13]:
optimal_weights = find_optimal_holdings(
    alpha_vector
)
print("Optimal weights: ", optimal_weights)
-----------------------------------------------------------------
           OSQP v0.6.2  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 3, constraints m = 0
          nnz(P) + nnz(A) = 0
settings: linear system solver = qdldl,
          eps_abs = 1.0e-03, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 500
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: on, time_limit: off

iter   objective    pri res    dua res    rho        time
   1  -3.9192e+06   0.00e+00   1.22e+00   1.00e-01   2.44e-04s
  25  -1.0000e+30   0.00e+00   1.22e+00   1.00e-01   6.17e-04s

status:               dual infeasible
number of iterations: 25
run time:             9.81e-04s
optimal rho estimate: 1.00e-06

Optimal weights:        0
0  None

As expected, without constraints, the problem is unbounded. We could achieve a smaller and smaller objective function value by taking larger and larger long positions on Stock C, and larger and larger short positions on Stock B. But this would increase risk and leverage. This is where our constraints come into play.

Create a Risk Model

We are going to create a PCA-based risk model for our 3 stocks. First, let's plot the returns data and their mean, just to have a sense for the shape of the data.

In [14]:
m = rets.mean()
In [15]:
# Trace for mean return vector m
mean_vec = np.vstack((np.full(3, 0), m.values)).T

hover_text2 = helper.generate_hover_text(mean_vec[0], mean_vec[1], mean_vec[2], 'Mean of Returns of Stock A', 'Mean of Returns of Stock B', 'Mean of Returns of Stock C', sig_figs = 7)

trace2 = go.Scatter3d(
    x=mean_vec[0],
    y=mean_vec[1],
    z=mean_vec[2],
    mode='lines+markers',
    marker=dict(
        size=4,
        color='red',
        opacity=0.9,
        symbol="diamond"
    ),
    name = 'mean daily return',
    text = hover_text2.flatten(),
    hoverinfo = 'text'
)


# Trace for data
hover_text3 = helper.generate_hover_text(rets['A'].values, rets['B'].values, rets['C'].values, 'Return of Stock A', 'Return of Stock B', 'Return of Stock C')

trace3 = go.Scatter3d(
    x=rets['A'].values,
    y=rets['B'].values,
    z=rets['C'].values,
    mode='markers',
    marker=dict(
        size=4,
        color='#7FB3D5',  
        opacity=0.3,
    ),
    name = 'daily returns',
    text = hover_text3.flatten(),
    hoverinfo = 'text'
)

data = [trace2, trace3]

layout = helper.create_standard_layout('Returns Data', 'Return of')
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig)

There are a few things to notice here.

  1. The vector of mean returns is very small compared to the volatility of the data! Admittedly, this is an artificial example with fake data, but this is also generally the case.
  2. The volatility of Stock C is larger than the volatility of Stock B, which is larger than the volatility of Stock A.

Given these data, we can already set up some expectations for what the results of PCA will look like. Remember that we expect the first PC to point in the direction of maximum variability of the data. Given the volatility of Stock C, we expect the first PC to point approximately parallel to the "Stock C axis". If we get results that are very different from this, we know that something is wrong with what we are doing. Let's now run the PCA algorithm, and see what we get.

First we'll mean-center the data.

In [16]:
rets = rets - m

Then we'll create a class that will fit our risk model.

In [17]:
from sklearn.decomposition import PCA

class RiskModelPCA():
    
    ANN_FACTOR = 252
    
    def __init__(self, num_factors):
        self._num_factors = num_factors
        self.num_stocks_ = None
        self.factor_betas_ = None
        self.factor_returns_ = None
        self.common_returns_ = None
        self.residuals_ = None
        self.factor_cov_matrix_ = None
        self.idio_var_matrix_ = None
        self.explained_variance_ratio_ = None

    def fit(self, returns):
        self.num_stocks_ = len(returns.columns)
        mod = PCA(n_components=self._num_factors, svd_solver='full')
        mod.fit(returns)
        
        self.factor_betas_ = pd.DataFrame(
            data=mod.components_.T,
            index=returns.columns
        )
        
        self.factor_returns_ = pd.DataFrame(
            data=mod.transform(rets),
            index=returns.index
        )
        
        self.explained_variance_ratio_ = mod.explained_variance_ratio_
        
        self.common_returns_ = pd.DataFrame(
            data=np.dot(self.factor_returns_, self.factor_betas_.T),
            index=returns.index
        )
        self.common_returns_.columns = returns.columns
        
        self.residuals_ = (returns - self.common_returns_)
        
        self.factor_cov_matrix_ = np.diag(
            self.factor_returns_.var(axis=0, ddof=1)*RiskModelPCA.ANN_FACTOR
        )
        
        self.idio_var_matrix_ = pd.DataFrame(
            data=np.diag(np.var(self.residuals_))*RiskModelPCA.ANN_FACTOR,
            index=returns.columns
        )
        
        self.idio_var_vector_ = pd.DataFrame(
            data=np.diag(self.idio_var_matrix_.values),
            index=returns.columns
        )
        
        self.idio_var_matrix_.columns = index=returns.columns

    def get_factor_exposures(self, weights):
        B = self.factor_betas_.loc[weights.index]
        return B.T.dot(weights)

Let's fit the risk model with 2 factors (i.e., we'll keep 2 PCs).

In [18]:
rm = RiskModelPCA(2) # create an instance of the class with 2 factors
rm.fit(rets) # fit the model on our data

Let's look at the factors (PCs).

In [19]:
rm.factor_betas_
Out[19]:
0 1
A -0.015530 0.059674
B -0.055588 -0.996725
C -0.998333 0.054570

The first PC points almost completely in the "Stock C" direction. The second points almost completely in the "Stock B" direction. This makes sense given the fact that we see most variability in the "Stock C" direction, and if we collapse that variability, we'd see the next most variability in the "Stock B" direction. Let's plot the PC vectors to make this more clear.

In [20]:
PC_scaler = 0.04 # The PC vectors have length 1, but this is larger than the scale of our data, so for visualization purposes, let's plot scaled-down versions of the PCs. 

# Trace for PC 0
pc0 = np.vstack((np.full(3, 0), rm.factor_betas_[0].values)).T*PC_scaler

hover_text4 = helper.generate_hover_text(pc0[0], pc0[1], pc0[2], 'Return of Stock A', 'Return of Stock B', 'Return of Stock C')

trace4 = go.Scatter3d(
    x=pc0[0],
    y=pc0[1],
    z=pc0[2],
    mode='lines+markers',
    marker=dict(
        size=4,
        color='#45B39D',
        opacity=0.9,
        symbol="diamond"
    ),
    line=dict(
        color='#45B39D',
        width=3
    ),
    name = 'PC 0',
    text = hover_text4.flatten(),
    hoverinfo = 'text'

)

# Trace for PC 1
pc1 = np.vstack((np.full(3, 0), rm.factor_betas_[1].values)).T*PC_scaler

hover_text5 = helper.generate_hover_text(pc1[0], pc1[1], pc1[2], 'Return of Stock A', 'Return of Stock B', 'Return of Stock C')

trace5 = go.Scatter3d(
    x=pc1[0],
    y=pc1[1],
    z=pc1[2],
    mode='lines+markers',
    marker=dict(
        size=4,
        color='#FFC300',
        opacity=0.9,
        symbol="diamond"
    ),
    line=dict(
        color='#FFC300',
        width=3
    ),
    name = 'PC 1',
    text = hover_text5.flatten(),
    hoverinfo = 'text'

)

# Trace for data
hover_text6 = helper.generate_hover_text(rets['A'].values, rets['B'].values, rets['C'].values, 'Return of Stock A', 'Return of Stock B', 'Return of Stock C')

trace6 = go.Scatter3d(
    x=rets['A'].values,
    y=rets['B'].values,
    z=rets['C'].values,
    mode='markers',
    marker=dict(
        size=4,
        color='#7FB3D5', 
        opacity=0.3,
    ),
    name = 'daily returns',
    text = hover_text6.flatten(),
    hoverinfo = 'text'
)

data = [trace4, trace5, trace6]

layout = helper.create_standard_layout('Returns Data with Factor (PC) Directions', 'Return of')
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig)

Let's look at the fraction of variance explained by the two factors we kept. This should confirm our intuition about the relative scale of variance explained by the first two factors.

In [21]:
plt.bar(np.arange(2), rm.explained_variance_ratio_, color=['#45B39D', '#FFC300']);
plt.title('% of Variance Explained by Each Factor');

Now let's look at the factor returns. We'll convert them to "price series" and visualize the evolution of these over time. Recall that these are the data expressed in the "factor" (PC) basis––the projections of the data vectors onto the factor dimensions.

In [22]:
(rm.factor_returns_ + 1).cumprod().plot(color=['#45B39D', '#FFC300'])
Out[22]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f1933b39d30>

Note that in this example, since the first factor is so close to the direction of the return on Stock C, the factor return (converted to a "price series") for factor 0 looks a lot like the price evolution of Stock C, but inverted. This is because the PCs are vectors that represent directions in space—they are only unique up to a sign change, and thus the projections onto them are unique up to a sign change.

In [23]:
prices.plot()
Out[23]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f1928081160>

Create Optimization Constraints

As we saw above, we need to constrain our optimization problem. Let's now create and plot the following constraints:

  • risk, based on our risk model
  • leverage
  • market neutral
  • factor exposure limits
  • individual weight limits

Risk

In [24]:
B = rm.factor_betas_
F = rm.factor_cov_matrix_
S = np.diag(rm.idio_var_vector_.values.flatten())

Using the $\mathbf{B}$, $\mathbf{F}$, and $\mathbf{S}$ matrices, let's write a function that takes in the portfolio weights and spits out the portfolio risk. Remember that the portfolio risk is calculated $\mathbf{x}^T(\mathbf{B}^T\mathbf{F}\mathbf{B} + \mathbf{S})\mathbf{x}$, but be careful to ensure the matrices are all oriented correctly.

In [26]:
def risk_fcn(x):
    """
    Calculate portfolio risk.

    Parameters
    ----------
    x : numpy array
        The vector of portfolio weights.

    Returns
    -------
    risk : float
        Portfolio risk.
    """
    #TODO: Implement function
    f = np.dot(B.values.T, x)
    risk = np.dot(f.T, np.dot(F, f)) + np.dot(x, np.dot(S, x))
    
    return risk
In [27]:
n_samples = 25
grid_max = 2.5
risk_grid, spacing, xv, yv, zv = helper.evaluate_fcn_on_grid(grid_max, n_samples, risk_fcn)

Now let's look at a plot of the value of the portfolio risk as a function of the portfolio weights in 3 dimensions.

In [28]:
hover_text = helper.generate_hover_text(xv, yv, zv, 'Weight on Stock A', 'Weight on Stock B', 'Weight on Stock C', fcn_values=risk_grid, fcn_label='Portfolio Risk')

trace7 = go.Scatter3d(
    x=xv.flatten(),
    y=yv.flatten(),
    z=zv.flatten(),
    mode='markers',
    marker=dict(
        size=4,
        color=risk_grid.flatten(), 
        colorscale='Viridis',
        opacity=0.3,
        showscale=True
    ),
    text = hover_text.flatten(),
    hoverinfo = 'text'
)

data = [trace7]
layout = helper.create_standard_layout('Portfolio Risk as Modeled by our PCA Risk Model', 'Weight on')
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig)

This plot uses a color at several points in 3-D space to represent the portfolio risk as a function of weight on each stock. Note that it increases with weight on each stock. This helps us visualize the risk function, but what does it look like to constrain risk? Let's visualize the boundary of the region within which risk is below a limiting value. This will help us visualize the shape of the space we search with optimization. We will do this for each of our constraints.

In [29]:
risk_data = helper.plot_iso_surface(risk_grid, 0.05, 25, 2.5, "Portfolio Risk = 0.05", '#F5B7B1', True)

Leverage

Now let's write a function that takes in the portfolio weights and spits out the portfolio leverage, $\sum_i|x_i|$. This will allow us to visualize the space that satisfies the constraint $\sum_i|x_i| \leq 1$.

In [30]:
def lev_fcn(x):
    """
    Calculate portfolio leverage.

    Parameters
    ----------
    x : numpy array
        The vector of portfolio weights.

    Returns
    -------
    leverage : float
        Portfolio risk.
    """
    #TODO: Implement function
    leverage = np.abs(x[0]) + np.abs(x[1]) + np.abs(x[2]) 
    
    return leverage

n_samples = 25
grid_max = 2.5

lev_grid, _, _, _, _ = helper.evaluate_fcn_on_grid(grid_max, n_samples, lev_fcn)
lev_data = helper.plot_iso_surface(lev_grid, 1.0, n_samples, grid_max, "Leverage Ratio = 1", '#1F618D', True)

Market Neutral

For the market neutral constraint, we constrain the sum of the weights. Let's write a function that calculates this quantity, given the portfolio weights. Then we can visualize the plane $\sum_i x_i = 0$.

In [31]:
def mn_fcn(x):
    """
    Calculate the sum of the portfolio weights.

    Parameters
    ----------
    x : numpy array
        The vector of portfolio weights.

    Returns
    -------
    sum_of_weights : float
        Sum of the portfolio weights.
    """
    #TODO: Implement function
    sum_of_weights = sum(x)
    
    return sum_of_weights

n_samples = 25
grid_max = 2.5

mn_grid, _, _, _, _ = helper.evaluate_fcn_on_grid(grid_max, n_samples, mn_fcn)
mn_data = helper.plot_iso_surface(mn_grid, 0, n_samples, grid_max, "Sum of Weights = 0", '#D35400', True)

Factor Exposures

Let's also calculate the factor exposures. Then we'll plot the planes defined by the factor exposure limits. The optimization will constrain the solution to lie between each pair of planes. We'll constrain each factor exposure to be between -0.1 and 0.1.

In [32]:
def fac_fcn(x):
    """
    Calculate portfolio factor exposures.

    Parameters
    ----------
    x : numpy array
        The vector of portfolio weights.

    Returns
    -------
    fac_exposures : numpy array
        Portfolio factor exposures.
    """
    #TODO: Implement function
    fac_exposures = np.dot(B.values.T, x)
    
    return fac_exposures

grid_max = 2.5
n_samples = 25

fac_grid, spacing, _, _, _ = helper.evaluate_fcn_on_grid(grid_max, n_samples, fac_fcn)

factor_limit_data = []

for factor in range(B.shape[1]):
    for l in range(2):
        mult_fac = l*2-1
        factor_limit_data.extend(helper.plot_iso_surface(fac_grid[:,:,:,factor], 0.1*mult_fac, n_samples, grid_max, 'Factor Exposure Limits', '#D2B4DE', False))

layout = helper.create_standard_layout('Factor Exposure Limits', 'Weight on')

fig = go.Figure(data=factor_limit_data, layout=layout)         
        
py.offline.iplot(fig)

Individual Weight Limits

Finally, let's look at the space inside the constraints on each individual weight. We'll constrain each individual weight to be between -0.55 and 0.55.

In [33]:
x_max = 0.55

x = np.array([-1, -1, 1, 1, -1, -1, 1, 1])*x_max
y = np.array([-1, 1, 1, -1, -1, 1, 1, -1])*x_max
z = np.array([-1, -1, -1, -1, 1, 1, 1, 1])*x_max
hover_text = helper.generate_hover_text(x, y, z, 'Weight on Stock A', 'Weight on Stock B', 'Weight on Stock C')


weight_data = [
    go.Mesh3d(
        x = x,
        y = y,
        z = z,
        colorscale = '#FCF3CF',
        intensity = np.full(8, 1),
        i = [7, 0, 0, 0, 4, 4, 6, 6, 4, 0, 3, 2],
        j = [3, 4, 1, 2, 5, 6, 5, 2, 0, 1, 6, 3],
        k = [0, 7, 2, 3, 6, 7, 1, 1, 5, 5, 7, 6],
        name='Weight on Stock B',
        showscale=False,
        opacity = 0.3,
        hoverinfo='none'
    )
]

trace = go.Scatter3d(
    x=x,
    y=y,
    z=z,
    mode='markers',
    marker=dict(
        size=6,
        opacity=0.0001,
        color='#BFB1A8', 
    ),
    text = hover_text.flatten(),
    hoverinfo = 'text',
    showlegend=False
)
    
weight_data = [weight_data[0], trace]    
    
layout = helper.create_standard_layout('Individual Weight Limits', 'Weight on')
fig = go.Figure(data=weight_data, layout=layout)
py.offline.iplot(fig)

Plot all the constraints on the same axes

Finally, let's visualize the space that satisfies all the constraints. It has to be inside all of the enclosed spaces, on the plane defined by the market neutral constraint, and between the planes defined by the factor exposure constraints.

In [34]:
data = risk_data
data.extend(lev_data)
data.extend(mn_data)
data.extend(factor_limit_data)
data.extend(weight_data)
layout = helper.create_standard_layout('Visualize Intersection of All Constraints', 'Weight on')
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig)

Optimization

Finally, let's define and run the optimization problem with the same objective function as above, but this time, with all the constraints we just discussed. Remember that the optimization is looking for the point in the space that satisfies all the constraints that minimizes the alpha function.

In [37]:
def find_optimal_holdings(
    alpha_vector,
    risk_model,
    risk_cap=0.05,
    factor_max=10.0,
    factor_min=-10.0,
    x_max=0.55,
    x_min=-0.55):
    
    """
    Define an optimization problem. It takes in several inputs and optimization parameters and outputs the
    optimized weights. The objective should minimize the objective -alpha*x, but also apply the risk, leverage,
    market neutral, factor exposure, and invidiual weight constraints.

    Parameters
    ----------
    alpha_vector : numpy array
        The alpha vector used in the optimization objective.
    risk_model : RiskModelPCA
        The risk model calculated above using PCA.
    risk_cap : float
        The limit to be placed on portfolio risk.
    factor_max : float
        The factor exposure upper limit.
    factor_min : float
        The factor exposure lower limit.
    x_max : float
        The individual weight upper limit.
    x_min : float
        The individual weight lower limit.

    Returns
    -------
    optimal_weights : numpy array
        The optimal portfolio weights.
    """
    #TODO: Implement function
    B = risk_model.factor_betas_
    F = risk_model.factor_cov_matrix_
    S = np.diag(risk_model.idio_var_vector_.values.flatten())
    
    x = cvx.Variable(len(alpha_vector))
    f = B.values.T*x
    
    risk = cvx.quad_form(f, F) + cvx.quad_form(x, S)
    
    obj = cvx.Minimize(-alpha_vector.values.flatten()*x)
    
    constraints = [
        sum(cvx.abs(x)) <= 1.0,
        sum(x) == 0,
        x <= x_max,
        x >= x_min,
        risk <= risk_cap,
        B.values.T*x <= factor_max,
        B.values.T*x >= factor_min
    ]
    
    prob = cvx.Problem(obj, constraints)
    prob.solve(verbose=True, max_iters=500)
    
    optimal_weights = np.asarray(x.value).flatten()
    
    return pd.DataFrame(data=optimal_weights, index=alpha_vector.index)
In [38]:
optimal_weights = find_optimal_holdings(
    alpha_vector,
    rm
)
In [39]:
optimal_weights
Out[39]:
0
A 3.857385e-11
B 5.000000e-01
C -5.000000e-01

The optimal weights are 0 weight on Stock A, 50% long on Stock B and 50% short on Stock C.

The solution for this notebook is here.

In [ ]: